#include "Distribution.h"

#include <limits>
#include <math.h>
#include <gsl/gsl_statistics.h>

#include <QDebug>

#define LOW_LIKELIHOOD 0.0000000000001

Distribution::Distribution() :
        randomNumberGenerator(gsl_rng_alloc(gsl_rng_ranlxd2))
{
}

Distribution* Distribution::sharedInstance()
{
    static Distribution distribution;
    return &distribution;
}

Distribution::~Distribution()
{
    gsl_rng_free(randomNumberGenerator);
}

double Distribution::randomUniform(double start, double end)
{
    return gsl_ran_flat(randomNumberGenerator, start, end);
}

double Distribution::randomNormal(double mean, double sigma)
{
    return mean + gsl_ran_gaussian(randomNumberGenerator, sigma);
}

double Distribution::randomGamma(double alpha, double beta)
{
    return gsl_ran_gamma(randomNumberGenerator, alpha, beta);
}

double Distribution::randomBernoulli(double p)
{
    return gsl_ran_bernoulli(randomNumberGenerator, p);
}

double Distribution::likelihoodNormal(double x, double mean, double sigma)
{
    return gsl_ran_gaussian_pdf(x - mean, sigma);
    //    double var = (x - mean) / sigma;
//    return 1.0/sqrt(2*M_PI) * exp(-0.5 * var * var) / sigma;
}

double Distribution::likelihoodGamma(double x, double alpha, double beta)
{
    double likelihood = gsl_ran_gamma_pdf(x, alpha, beta);
    if (likelihood <= LOW_LIKELIHOOD)
        likelihood = LOW_LIKELIHOOD;
    return likelihood;
}

double Distribution::likelihoodInverseGamma(double x, double alpha, double beta)
{
    double result = alpha * log(beta);
    result -= ((alpha + 1) * x * log(x) + x * log(this->gamma(alpha)) + beta)/x;
    result = exp(result);
    return result;
}

double Distribution::likelihoodBernoulli(int x, double p)
{
    double likelihood = gsl_ran_bernoulli_pdf(x, p);
    if (likelihood <= LOW_LIKELIHOOD)
        likelihood = LOW_LIKELIHOOD;
    return likelihood;
}

double Distribution::likelihoodBeta(double x, double alpha, double beta)
{
    double likelihood = gsl_ran_beta_pdf(x, alpha, beta);
    if (likelihood <= LOW_LIKELIHOOD)
        likelihood = LOW_LIKELIHOOD;
    return likelihood;
}

double Distribution::likelihoodPoisson(int x, double rate)
{
    double likelihood = gsl_ran_poisson_pdf(x, rate);
    if (likelihood <= LOW_LIKELIHOOD)
        likelihood = LOW_LIKELIHOOD;

    return likelihood;
}

double Distribution::mean(QVector<double> &data)
{
    return gsl_stats_mean(data.data(), 1, data.count());
}

double Distribution::variance(QVector<double> &data)
{
    return gsl_stats_variance(data.data(), 1, data.count());
}

double Distribution::minimum(QVector<double>& data)
{
    double min = std::numeric_limits<double>::max();
    foreach (double each, data) {
        if (each < min) {
            min = each;
        }
    }

    return min;
}

double Distribution::maximum(QVector<double>& data)
{
    double max = std::numeric_limits<double>::min();
    foreach (double each, data) {
        if (each > max) {
            max = each;
        }
    }

    return max;
}

QVector<QPointF> Distribution::histogramLine(QVector<double>& data, int bucketCount)
{
    double min = this->minimum(data);
    double max = this->maximum(data);

    double stride = fabs(max - min) / (double) (bucketCount);
    QVector<int> histogram = this->histogram(data, bucketCount);
    QVector<QPointF> histLine;
    double x = min;
    for (int y = 0; y < histogram.count(); y++, x+=stride) {
        QPointF point(x, histogram[y]);
        histLine.append(point);
    }

    return histLine;
}

QVector<int> Distribution::histogram(QVector<double> &data, int bucketCount)
{
    double min = this->minimum(data);
    double max = this->maximum(data);

    double stride = fabs(max - min) / (double) (bucketCount);
    QVector<int> histogram;
    histogram.fill(0, bucketCount + 1);
    foreach (double each, data) {
        double index = each/stride - min/stride;
//        if (index >= histogram.count())
//            qDebug() << index;
        histogram[(int) index]++;
    }

    return histogram;
}


double Distribution::gamma(double x)
{
    //  gamma.cpp -- computation of gamma function.
    //      Algorithms and coefficient values from "Computation of Special
    //      Functions", Zhang and Jin, John Wiley and Sons, 1996.
    //
    //  (C) 2003, C. Bond. All rights reserved.
    //
    // Returns gamma function of argument 'x'.
    //
    // NOTE: Returns 1e308 if argument is a negative integer or 0,
    //      or if argument exceeds 171.
    //
    int i,k,m;
    double ga,gr,r,z;

    static double g[] = {
        1.0,
        0.5772156649015329,
       -0.6558780715202538,
       -0.420026350340952e-1,
        0.1665386113822915,
       -0.421977345555443e-1,
       -0.9621971527877e-2,
        0.7218943246663e-2,
       -0.11651675918591e-2,
       -0.2152416741149e-3,
        0.1280502823882e-3,
       -0.201348547807e-4,
       -0.12504934821e-5,
        0.1133027232e-5,
       -0.2056338417e-6,
        0.6116095e-8,
        0.50020075e-8,
       -0.11812746e-8,
        0.1043427e-9,
        0.77823e-11,
       -0.36968e-11,
        0.51e-12,
       -0.206e-13,
       -0.54e-14,
        0.14e-14};

    if (x > 171.0) return 1e308;    // This value is an overflow flag.
    if (x == (int)x) {
        if (x > 0.0) {
            ga = 1.0;               // use factorial
            for (i=2;i<x;i++) {
               ga *= i;
            }
         }
         else
            ga = 1e308;
     }
     else {
        if (fabs(x) > 1.0) {
            z = fabs(x);
            m = (int)z;
            r = 1.0;
            for (k=1;k<=m;k++) {
                r *= (z-k);
            }
            z -= m;
        }
        else
            z = x;
        gr = g[24];
        for (k=23;k>=0;k--) {
            gr = gr*z+g[k];
        }
        ga = 1.0/(gr*z);
        if (fabs(x) > 1.0) {
            ga *= r;
            if (x < 0.0) {
                ga = -M_PI/(x*ga*sin(M_PI*x));
            }
        }
    }
    return ga;
}
